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ABSTRACT 

We present a new method for modeling inhomogeneous cosmic reionization on large scales. Uti- 
lizing high-resolution radiation-hydrodynamic simulations with 2048 3 dark matter particles, 2048 3 
gas cells, and 17 billion adaptive rays in a L = 100 Mpc/h box, we show that the density and 
rcionization-rcdshift fields are highly correlated on large scales (> 1 Mpc/h). This correlation can 
be statistically represented by a scale-dependent linear bias. We construct a parametric function for 
the bias, which is then used to filter any large-scale density field to derive the corresponding spatially 
varying rcionization-rcdshift field. The parametric model has three free parameters which can be 
reduced to one free parameter when we fit the two bias parameters to simulations results. We can 
differentiate degenerate combinations of the bias parameters by combining results for the global ion- 
ization histories and correlation length between ionized regions. Unlike previous semi-analytic models, 
the evolution of the reionization-redshift field in our model is directly compared cell by cell against 
simulations and preforms well in all tests. Our model maps the high-resolution, intermediate-volume 
radiation-hydrodynamic simulations onto lower-resolution, larger-volume N-body simulations (> 2 
Gpc/h) in order to make mock observations and theoretical predictions. 

Subject headings: Cosmology: Theory — Galaxies: Clusters: General — Large-Scale Structure of 
Universe — Methods: Numerical 



1. INTRODUCTION 

When the first stars and galaxies turned on and began 
ionizing the surrounding cold and neutral hydrogen of the 
intcrgalatic medium (IGM), they started the phase tran- 
sition of the Universe known as th e Epoch of Reionization 
fEoR: lLoeb fe Furlanettol 120131) . This inhomo geneous 
process of reionization leaves two, among others, possible 
observable sources, the neutral hydrogen atoms and the 
ionized electrons. It is necessary that precise theoretical 
models of EoR on Gpc scales are used to interpret the 
information from the observable imprints left by these 
sources in order to gain insight into and understanding 
of th e first stars and the initial stages of galaxy forma tion 
(e.g. lFurlanetto et al.|[200llMorales fc Wvithe|[2010l and 
references therein). 

Neutral hydrogen atoms are observed in both absorp- 
tion and emission. Current constraints from absorp- 
tion measurements come from observations zero trans- 
mission of rest-frame Lya flux at z > 6 in spectra of 
high redshift qu asars, which su ggest that EoR is com- 
pleted by z ~ 6 (jFan et al.H2006f ). although it is possible 
for these constraints to be consistent with reionization 
completing at a high er redshift (e.g. IQh fc Furlanettol 
[2005l:lLidz et~aL|[2006h . The neutral hydro gen emission is 
observable through the redshifted 21cm si gnal that orig- 



inates from its hyperfine transition (e.g. Scott fc Reesl 
[ll)90t ISha7er et al.lfl999t IZaldarriaga et al.ll2004l) . There 



are several experiments currently searching for the 21cm 
signal at z > 6, such as, the Murchison Wide Field 
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Array (MWA 4 ; IBowman et "all [20051), the Giant Meter- 
wave Telescope (GMRT 5 : |Pen et al.H2009l), th e Low Fre- 
quency Array (LOFAR 6 : lHarker et all 120101 ), and The 
Precision Array for Probing th e Epoch of Reionization 
fPAPER 7 : iParsons et all 120101) . Using the 21cm signal, 
the experiment EDGES 8 reported a lower limit to du- 
ratio n of reionization, A z > 0.06 (IBowman fc Rogers! 
I2012I ). Future 21 cm experiments like the proposed 
MeeRKAT 9 (|Booth et al]|2009|) and the Square Kilome- 
ter Array fSKA 10 : [Mellema et all I2012D both have the 
potential to measure this signal across several frequen- 
cies and provide tomographic information on EoR. 

The free electrons are observed through their scat- 
tering of cosmic microwave background (CMB) pho- 
tons. This scattering can be seen on large scales in 
polarization CMB or on small scales in the CMB sec- 
ondary anisotropies, such as kin etic Sunyaev Zel ' dovich 
(kSZ) sign a l fro m the EoR (iGruzinov fc Hv3 [l998l: 



2003; 



2007 



Zahn et al. 



Knox et all Il998t IValageas et all 120011: iSantos et al 



2005t IMcQuinn et all 120051: llliev et a! 



Mesinger et alJ 120111 ). Polarization measurements 



of th e CMB place a con straint on the optical depth to 
EoR |L arson et "all 120 111 ). Assuming a step function or 
a hyperbolic tangent function for the ionization history, 
the 7-year WMAP data implies that the reionization- 
redshift is 10.5 ± 1.2 (68% CL). Constraints have also 
come from multi-frequency high resolution CMB exper- 
iments that measure the power spectrum of CMB sec- 
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Fig. 1. — Slices from our high resolution RadHydro simulation for a model of reionization that occurs "late" with a midpoint of z = 8 
and is finished by z fa 6.9. The dimensions are 100 Mpc/h x 100 Mpc/h with a thickness of ~100 kpc//i comoving. Left: The density 
field, p(x)/p. Right: The reionization-redshift field, zre(x). Large-scale, overdense regions near sources are generally ionized earlier than 
large-scale, undcrdense regions far from sources. 



ondary anisotropics to great precision, where contribu- 
tions to kSZ power fro m EOR are the lar gest. The South 
Pole Telescope (SPT n : IZahn et al.ll20T2l) placed a model 
dependent upper limit on the duration of reionization 
from their multifrcquency measurements of the high I 
power spectrum, and future results from the Atacama 
Cosmology Telescope (ACT 12 ) are expected to place sim- 
ilar constraints. The next generation high resolution 
CMB experiments ACT with polarization (ACT-pol) and 
South Pole Telescope with polarization (SPT-pol) will 
precisely measurement the secondary anisotropics of the 
CMB in both temperature and polarization, which will 
provide tighter constraints on EoR. 

For the EoR experiments listed above and future ones, 
the amount of understanding gained on these first ioniz- 
ing sources and the initial stages of galaxy evolution will 
depend upon the accuracy of the theoretical models for 
EoR. The main challenge in EoR theory is providing an 
accurate model of the IGM, the sources and the sinks of 
ionizing photons, while having the a large enough vol- 
ume > l(Gpc//i) 3 to statistically sample the HI regions 
and construct mock observations on the angular scales 
required by the current and future EoR experiments. 

There are two standard approaches to model EoR, 
radiative transfer simulations with various imple- 
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a region is fully ionized if the simple relation, 
coil > 1 is satisfied. Here £ is an efficiency parameter 
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els 

and -Fcoii is the collapse fractio n, which is calcul ated 
via the excursion set formalism (|Bond et al.l 119911) , or 
appli ed to three dimen sional realization of a density field 
(e.g. IZahn et al.l [20051 ). Semi-analytic models capture 
the generic properties of EoR, but in order to capture 
the complex non-linear, and non-Gaussian nature of 
EoR radiative transfer simulations are required. 

The advantages of the current full hydrodynamic, high 
resolution simulations with radiative transfer (imple- 
mented either in post processing or during the sim- 
ulation) is that they probe the relevant scales to re- 
solve sources of ionizing photons and their sinks, then 
trace these photons through an inhomogeneous IGM 
(|Trac fc G nedin 2011)). However, full hydrodynamic sim- 
ulations with radiative transfer on large enough scales 
to capture a representative sample of ionizing sources 
and with enough small scale resolution to also capture 
all the physics of reionization are currently not possi- 
ble due to the overwhelming computational demands of 
such calculations. Thus, all of the simulations to date 
ha ve been rest r icted to smaller box-sizes. Recent work 
by IZahn et al.l (|2011l) ran several convergence tests be- 
tween these two types of EoR models. For all the models 
in their study, they found that the results from the mod- 
els are within tens of percent of each other. Although in 
these comparisons the parameters of semi-analytic mod- 
els were adjusted to match the ionization fractions of the 
simulations at the redshifts of interest. 

In this paper, we present a substantially more accu- 
rate semi-analytical model that is statistically informed 
by simulations with radiative transfer and hydrodynam- 
ics. The implementation of this model is fast, versatile 



A Parametric model for Cosmic Reionization 



3 




Fig. 2. — The correlation, f'mz> and bias, b mz (fc), calculated from simulations as a function of scale (A:). Left: Results for r mz from two 
different simulated reionization scenarios are shown by the red line(early reionization) and blue line(late reionization) with a resolution of 
L = 100 Mpc//i box length within 2048 cells (solid lines). In both scenarios r mz decreases rapidly on scales below 1 Mpc//i and the high 
resolution simulations show higher correlations to smaller scales. The overall shape of r mz is largely independent of whether reionization 
occurred early or late and is completely independent of the ionization threshold chosen (90 - 50%) for the zre(x) field (solid blue and 
dashed orange lines, respectively). Right: Results for fc mz (fc) from simulations (solid lines), the fiducial model, which is the best fit k a and 
a parameters to fc mz (fc) from simulations (dotted line), and the b mz (fc) functions that represent the long and short duration reionization is 
the orange and light green dashed line, respectively. The horizontal grey line represents the analytical prediction from spherical collapse. 



and easily applied to large N-body simulations, thus it 
can be scaled up to the large volumes required by the 
current and future EoR experiments without loss of ac- 
curacy. This paper is the first in a series that explore 
EoR observables produced via our model. We focu s 
on CMB relate d obse rvables in iNataraian et al.l (12012D: 
iBattaelia et al.1 (|2012t ) and the 21cm in lLa Plante et al.1 
(|2012[ ). In section [21 we present our fast semi-analytical 
model and the simulations it is calibrated on. Section [3] 
compares the model to the simulations on a cell by cell 
bases. We show results on the global reionization his- 
tory and the typical correlation between ionized regions 
in section |U We compare to previous work and discuss 
caveats to our model in section [S] and conclude in sec- 
tion [5] Throughout the paper, we adopt the concor- 
dance cosmological parameters: O m = 0.27, fl\ = 0.73, 
Q b = 0.045, h = 0.7, n s = 0.96, and cr 8 = 0.80. 

2. METHODOLOGY 

We present a novel semi-analytic model for calculating 
the evolution of the 3-dimensional ionization field in large 
volumes > (Gpc//i) 3 , which is currently not attainable 
with direct simulations. Using radiation-hydrodynamic 
simulations, we demonstrate that the rcdshift at which 
a volume-element is ionized can be calculated by fil- 
tering a nonlinear density field with a simple paramet- 
ric function. Our method can be used to map high- 
resolution, intermediate- volume radiation-hydrodynamic 
simulations onto lower-resolution, larger-volume N-body 
simulations in order to make mock observations and the- 
oretical predictions. In addition, the model parameters 
can be varied away from the fiducial values in order to 
explore the reionization parameter space (e.g. the timing 
and duration of the EoR) . 

2.1. Hydrodynamic Simulations with Radiative 
Transfer 



We adopt the hybrid approach in simulating cosmi c 
reionization previously described in iTrac et all (|2008l) . 
First, a high-resolution N-body simulation is used to 
evolve the matter distribution and track the formation of 
dark matter halos. The resulting halo catalogs are used 
to develop a subgrid model for high-rcdshift radiation 
sources. Second, direct RadHydro (radiative transfer + 
hydrodynamic + N-body) simulations are used to simul- 
taneously solve the coupled evolution of the dark matter, 
baryons, and radiation. 

A particle-particlc-particle-mesh (P 3 M) code is used to 
run a high-resolution simulation with 3072 3 dark matter 
particles in a 100 Mpc/h comoving box. A spherical over- 
density halo finder is used on the fly to identify collapsed 
dark matter halos with average densities equal to 200 
times the average cosmic density. With a particle mass 
resolution of 2.58 x 10 6 M & /h, we can reliably locate dark 
matter halos down to the atomic cooling limit (T ~ 10 4 
K, M ~ 10 8 M & /h). The halos are populated with ra- 
diation sources and the ionizing photon production rate 
N*,(M,z) is calcula ted using the halo model described in 
ITrac fc Cenl (I2007T) . 

The RadHydro code combines a cosmological hydro- 
dynamic code (moving frame hyd rodynamics + particle- 
mesh N-body; ITrac fc Pen! I2004T ) wit h an adaptive ray - 
tracing radiative transfer algorithm ([Trac fc Cenl 120071 ). 
The raytracing algorithm has adaptive splitting and 
merging and utilizes a two-level radiative transfer (RT) 
grid scheme to obtain better resolution and scaling. We 
have run two moderate-resolution RadHydro simulations 
each with 2048 3 dark matter particles, 2048 3 gas cells, 
and 17 billion adaptive rays. The RadHydro and N-body 
simulations were run using the Blacklight supercomputer 
at the Pittsburgh Supercomputing Center (PSC). 

In the first RadHydro simulation, reionization occurs 
earlier with a midpoint of z w 10 (mass and volume- 
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Fig. 3. — A comparison of the zre(x) fields (top panels) and the ionization fraction field at z = 8.1 (bottom panels) and between our 
RadHydro simulation with a late reionization scenario (left panels) and model (right panels). The color scales in the top panels illustrate 
the redshift at which these cells reionize and the black regions in the lower panels correspond to non-ionized regions, while the white regions 
are ionized z = 8.1. The simulations were degraded down to a resolution of 1 Mpc//i (see Sec. [3] for averaging details) and the slices 
thicknesses arel Mpc/h. In both upper and lower panels the model captures the same structures as the RadHydro simulation on large 
scales. The zre( x ) fields differ on small scales in the void and filament regions, where the our method tends to predict earlier redshifts of 
reionization. 



weighted ionization fractions are 0.5) and is effec- 
tively completed by z « 8.7 (radiation filling factor of 
the radiation- hydrodynamic grid reaches unity). This 
early reionization model has a Thomson optical depth 
for electron scattering r « 0.088, which is in good agree- 
ment with current observational constraints. From the 
WMAP 7-year results, the Thomson optical depth is 
r = 0.088 ± 0.015 assuming instantaneous reionization 
and r = 0.087 ± 0.015 if the wid th of reionization is 
allowed to vary ([Larson et al.ll20ilT ) . In the second simu- 
lation, reionization occurs later with a midpoint of z = 8 
and finished by z ss 6.9. The late model has r w 0.067, 
which is lower but within 2a of the WMAP best-fit value. 



In the two basic simulations, the radiative trans- 
fer of the ionizing photons proceeded such that large- 
scale, overdense regions near sources are generally ion- 
ized earlier than large-scale, underdense regions far from 
sources. The subgrid model for sources only included 
high- redshift galaxies with Population II stars (jSchaererl 
l200l since they are expected to provide the dominant 
contribution to the ionizing photon budget. The simula- 
tions do not include reionization by Population III stars 
or X-ray sources, which primarily affect th e earliest phase 
of the EoR (e.g. iFurlanetto et all l2006f ) . We also ne- 
glect additional clumping and self-shieldin g of small-scale 
dense absorbers such as mini-halos (e.g. IShapiro et all 
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Fig. 4. — A comparison of the cumulative distribution function 
(CDF) for the zre(x) field between RadHydro simulations (early 
and late reionization scenarios, represented by the red and blue 
lines, respectively) and our model predictions constructed from 
the same density fields as the simulations. Our models deviate 
from the RadHydro simulations at the beginning and end of reion- 
ization, since the models are constructed for densities fields at z. 
Overall this accounts for percent level changes in the duration of 
rcionization, the exact amount will depend on the definition for the 
duration. 



120041 ) or Lyman limit systems (e.g. lGnedin fc Fanll2006l ). 
We will explore other reionization scenarios using differ- 
ent models for sources and sinks in future work. 

2.2. Semi- analytic model 

For every cell in the RadHydro simulations, we 
record the rcdshift at which it ionizes and construct a 
rcionization-redshift field, zre(x). Figure Q] shows a slice 
of the the density and the reionization-redshift fields for 
the late reionization scenario. For the purposes of this 
work, regions that are greater than 90% ionized are con- 
sidered to be finished reionizing, but we obtain nearly 
identical statistical quantities when we chose 50% (cf. 
Fig [2]) . Defining a threshold for when cells are consid- 
ered to be ionized is arbitrary, since there is a sharp 
transition between the onset of ionization and comple- 
tion within a cell. Most cells approach 100% ionization 
but never reach it due to recombinat ions. Our results 
are consistent with lTrac fc Cenl (|2007t ) who used a 99% 
ionization threshold. 

Figure Q] shows the over-density field at z = 8 and 
the corresponding rcionization field from the simula- 
tions, which clearly illustrates that the reionization- 
redshift is associated with the density. The two fields 
arc highly correlated since large-scale over-densities near 
sources are generally ionized earlier than large-scale 
unde r-dense regions far from sources (|Barkana &: Loebl 
12004) . Other semi-analytical approaches implicitly in- 
voke this association when constructing th eir models 
(e.g. iFurlanetto eFall 121)041 : iZahn et al.1 [20051 ) . and Fig- 
ure Q] illustrates that this assumption is fairly accurate 
down to Mpc scales. Our method quantifies the corre- 
lation between the reionization-redshift and density in 
our RadHydro simulations and uses this correlation and 
its bias to construct reionization fields from any density 
field. First, we define the following fluctuations fields 
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Fig. 5. — The 2D probability distribution function (PDF) be- 
tween the zre(x) fields from the RadHydro simulation for late 
reionization and our model predictions constructed from the same 
density field compared on the cell by cell level. The contours repre- 
sent the cumulative probability distribution function for 50% and 
90% of the cells in this comparison. Our model preforms extremely 
well around z = z and slightly biases the zre (x) field to larger red- 
shift values. 



^m(x) 



p(x) - p 



and 



4(x) 



[1 + zre(x)] - [1 + z) 
1 + z 



(1) 
(2) 



where p is the mean matter density and z is the mean 
value for the zre(x) field. This construction has the 
advantage that we removed mean redshift dependences 
of zre(x) and these fields arc dimcnsionless. Then 
we quantify their correlation using two point statistics 
(SiSj) = Pij{k) in Fourier space (i.e. power spectra and 
cross spectra) assuming isotropy and we calculate their 
cross correlation 



Pita {k) 



A/Pzz(fc)Pmm(fc) 



and linear bias, 



1 P™(k) 

Pmm{k) 



(3) 



(4) 



Figure [2] show that these fields are highly correlated on 
scales above 1 Mpc/h, and on scales below 1 Mpc/h the 
correlation decreases. Similar cross correlati ons between 
ioniza tion fields and density are calculated in lZahn et al.l 
(| 2 1 If ) for fixed redshifts. As a result of r mz (fc) being 
highly correlated on scales above 1 Mpc//i, just knowing 
the bias between these two fields allows one to construct 
cither field from the other by filtering the initial field 
with the bias. The newly constructed field will statisti- 
cally match the results from the simulations. We chose to 
calculate a the linear bias (cf. Eq. [4} , since this bias does 
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well to match the zre(x) from the simulations on scales 
>1 Mpc/h (cf. Fig. [3]). Going to higher order correla- 
tions and bias models should improve this comparison, 
but is not necessary to capture the global features at the 
accuracy we are interested in; such extensions are left 
for future work. We find that the linear bias can be fit 
by a three parameter function (cf. Fig. [5]). The simple 
functional form for the bias factor is, 



b mz (k) 



(1 + k/k ) c 



(5) 



with the three free parameters being, b Q , k oy and a. A 
simple description for the parameters is as follows: b a 
is the bias amplitude, k a is the scale threshold, and a 
is the asymptotic exponent. We fit k a , and a to results 
from the simulations using least squares fitting with the 
correlation function (r mz ) as weight (cf. Fig. [2]). This 
weighting emphasizes the scales where S z and 5 m are 
highly correlated and down weights the small scales. The 
best fit values for the b mz (k) are k a = 0.185 Mpc/h and 
a = 0.564. Hereafter, we refer to these parameter val- 
ues as the fiducial bias parameters. We show in Figure [2] 
that our parametric function with these fiducial param- 
eters is comparable to calculated simulation bias. Also 
we explore the parameter space of k Q and a and allow 
them to vary about the fiducial parameters. 

The third parameter, b Q , which is the bias amplitude 
on the largest scales is not fit for, since the simulations 
that we ran have finite box sizes and these scales ex- 
tend beyond their box size. Figure [2] shows that b mz (k) 
is asymptotically approaching some value for 6 , but fit- 
ting for this value would be inaccurate and degenerate 
with k Q and a. Instead we ref er to previous analytic 
work by [Barkana fc Loebl (|2004f ) where this large-scale 
bias is derived using the excursion set formalism (i.e . 
extended press-schechter formalism; iBond et al.lll99ll ). 
They show that on these large-scales the differences in 
the rcdshift of collapse and rcioinization for various over- 
densites arc related via, 



S z {k^0) = 



S m (k -> 0) 



(6) 



Here S c = 1.68 is the critical over-density threshold. 
Throughout this work we use the value 1/S C for b Q = 
0.593. 

2.3. Implementation 
We used a particle-particle-particle-mesh (P 3 M) N- 
body code to evolve 2048 3 dark matter particles in a 2 
Gpc/h box to generate the over-density fields, <5 m , down 
to z = 5.5. We convolve the S m with a filter consist- 
ing three elements: (1) A cubical top hat filter, S(fc), 
which deconvolves the smoothing used to construct 5 m 
from simulation, (2) A Fourier transform of a real space 
top hat filter 0(fc), which smoothes <5 m to resolution of 1 
Mpc/h, and (3) The bias function from Equation[5] The 
assembled filter takes this form 



W z {k) 



b mz (k)e(k) 



(7) 




and we apply this filter at z. We Fourier transform back 
to real space and convert the newly constructed S z field 



Fig. 6. — The ionization fraction as a function of redshift, x e (z), 
for three fiducial parameter models with 2 = 8, 10, 12 (green, red, 
and blue lines, respectively) and the long and short duration reion- 
ization models (orange and light green dashed line, respectively). 
The fiducial models have similar ionization histories they are just 
shifted according to z. 

to zre(x) field using Eq.[5]with the same z as the density 
field, here z essentially sets the midpoint of reionization. 
Thus, we have a complete ionization history for the den- 
sity field used and from 2re(x) we can construct maps 
or 3D realization at a given redshift. The entire process 
is quick, since it only requires Fourier transforms, which 
can be done with any fast Fourier transform software and 
scale like NLog(N), where N is the number of cells. We 
found that I/O of these simulation data files dominates 
the total time when producing a zre(x) field. This al- 
lows one to easily explore the parameter space of fc Q , a, 
and z. 

3. COMPARISON TO SIMULATIONS 

We compare the results from the RadHydro simula- 
tions for the reionization- rcdshift fields, zre(x), to our 
model predictions that were constructed fie ld from the 
same simulated density fields. Similar to iZahn et al.1 
(|2011l) . we also compare their ionization fields at fixed 
redshifts. Unlike previous comparisons between simula- 
tions and semi-analytic models, our comparisons of the 
•zre(x) fields are computed at the cell by cell level and 
our model is not tuned to match a particular ionization 
fraction at a given redshift. Our comparisons go beyond 
bubble size and distribution tests and looks at the red- 
shift evolution of these quantities as well. In order to 
make these direct comparisons, the simulations must be 
smoothed to the scales on which we apply our model 
> lMpc//i, since they are at a higher resolution than 
1 Mpc//i. The simulations are smoothed by convolving 
them with a cubical tophat filter, which degrades the res- 
olution of their density and zre(x) fields to ~ 1 Mpc/h. 
After smoothing the simulations we find that the same 
structures in the zre (x) field slice seen in Figure [T] are 
still visible in Figure El 

Our model predictions for zre(x) are constructed on 
the density fields that are smoothed in the same way as 
the zre(x) fields, but the 6 mz (fc) used is calculated di- 
rectly from the original simulation. The top panels of 
Figure [3] show that the zre(x) fields from the Simula- 
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Fig. 7. — Parameter grids of k and a for the width of reionization, A z . In the left panel A z = z(x c = 25%) — z(x c = 75%) and in the 
right panel A z = z(x e = 5%) — z(x e = 95%). Increasing k a extends the duration of rcionization, while increasing a shortens the duration 
of reionization. The trends for A z (fc , a) are independent of the definition of A z . There are degenerate combinations of k and a that give 
the same values for A z . 



tions and our model agree on large scales in their spa- 
tial distribution, structure and evolution. Additionally, 
this illustrates that on smaller scales there are slight dis- 
agreements in the filament and void regions as well as 
the shape of the structures produced from the models, 
which are more filamentary than simulations. We find 
that the ionization fields model at z =8.1, which approxi- 
mately corresponds to 50% ionization in the zre (x) field, 
reproduces the simulations results at the same redshift 
(cf. Fig. [3]). The difference found between them are at- 
tributed to slightly different mean ionization fractions at 
z =8.1 and the steep changes in these fractions around 
this redshift. All the differences listed are attributed to 
the fact that the correlation is not exactly one at these 
scales, so using a model with a scale-dependent linear 
bias will not be perfect. 

For a quantitative comparison of our model we first 
compare simulation and model results for the cumula- 
tive distribution functions (CDF) of zre(x) for an early 
and late scenario of reionization. Figure |4] shows that our 
model traces the simulations results well and only devi- 
ates from simulations at the very beginning and end of 
rcionization, although this deviation is small. This cor- 
responds to our model having inaccuracies in the densest 
regions and in the voids, which is illustrated previously 
in Figured The zre(x) values where the CDF(z) = 0.5 
in our models show small offsets compared to the sim- 
ulations zre(x) values for the same point in the CDF. 
We preform a cell by cell comparison of the zre (x) fields 
between simulations and our model predictions, which 
are the most stringent test. In Figure [5] we show the 2D 
probability distribution function (by the blue color scale) 
and cumulative distribution function (red contours) for 
the late reionization scenario. A majority of the cells fall 
along the one to one line in this comparison with some 
scatter and they are clustered around z of this simulation. 
We find that 90% of the cells from the model arc within ~ 



10% of the simulation value for their 2re(x) values. Af- 
ter preforming these stringent cell by cell test we did not 
find any large systematic biases in the mean values or the 
variations, but these small differences mentioned above 
are negligible compared to the uncertainties in physics of 
reionization. We emphasize that no two radiative trans- 
fer hydrodynamic simulations nor semi-analytical models 
have ever been compared by such stringent tests. Thus, 
our method of filtering the density field with a scale- 
dependent linear bias is a sufficient description of reion- 
ization especially on scales larger than 1 Mpc/h. 

4. RESULTS 

In this section, we present results on two physical as- 
pects of the model, the ionization history and the correla- 
tion length between ionized regions, which is a proxy for 
bubble size. We explore the parameter space of Eq. [5] for 
aspects, which provides physical understanding of zre(x) 
field in relation to the parameters k a , a, and z. 

4.1. Ionization History 

We show the results for the mass weighted ionization 
history, x e , from our semi- analytic model for various pa- 
rameter values of k Q , ex, and z in Fig. [6] In this paper we 
present only the results for the mass weighted x e , since 
both the mass and volume weighted x c have comparable 
redshift evolution with the slight difference being that the 
mass weighted x e increases faster at higher redshifts than 
the volume weighted x e . The parameter z approximately 
sets the redshift where x c = 50%. We find that for a 
fixed k and a the shape of x e (z) about z is essentially 
independent of the value z and the percent differences 
between z = 8, 10, 12 are below 10% for x c > 10% (cf. 
Fig[S]). The redshift evolution of the ionization history is 
set by k a and a, where increasing a shortens the duration 
of reionization and increasing k a lengthens the duration 
of reionization (cf. Fig. [7J. When we increase a, we ex- 
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Fig. 8. — Left: The 6 mz (fc) functions for two examples of degenerate pair combinations of k a and a with similar A z at fixed z = 10 
(cf. Fig[7)l. Right: The x e (z) corresponding to the same pair combinations. Even though the 6 mz (fc) functions of these degenerate pair 
combinations are different, they have similar x B (z). 



tend 6 zm to small scales. Thus, we increase the variance 
in S z and the duration of reionization, while decreasing 
k has opposite effect. These effects can magnified by 
increasing one parameter while decreasing the other or 
vice versa. 

Similar to previous work (e.g. iZahn et al.l I2012H . we 
define two measures for the duration of reionization, 
A z = z(x c = 25%) — z(x c = 75%) and A z = z(x c = 
5%) — z(x c — 95%). In general, semi-analytic models and 
simulations exclude the early and late times of reioniza- 
tion when defining A z , since it is difficult to capture the 
small scale physical processes at these times. Figure [7] 
shows how k Q and a affect A z . The trends in x from 
varying the values k a and a are the same for A z and 
independent of the definition for A z . Currently there 
are upper limits for A z CMB small sca le measurements 
(|Zahn et al.ll2012t iMesinger et al.ll2012l ). independent of 
definition and lower limits for A y fr om global 21 cm ob- 
servations (jBowman fc Rogers 2012). With a few excep- 
tions, all our models fall well within the constrained re- 
gion of parameter space for A z from these observations, 
although, we emphasize that the conversion from obser- 
vation to A z is model dependent. 

There are degeneracies in both x c (z) and A z in our 
parametric bias model. Figure |8] shows two examples 
of parameters pairs values of a = 0.2, k a = 0.1 and 
a = 0.6, k Q = 0.7 that have different bias functions, but 
they have similar ionization histories (cf. Fig. [5]) and A z 
values (cf. Fig. [7J) . These degeneracies are a problem if 
one wants to relate any x e (z) observable to the underly- 
ing parameters of our model. It is possible to use higher 
order statistics than A z , i.e. beyond the variance, to 
differentiate between these degenerate models, but this 
does not add physical understanding of how a and k 
affect EoR. Any measures or proxies for the typical sizes 
of ionizing regions will differentiate between the degen- 
erate pairs of a and k a , while providing a more physical 
understanding of the impact these parameters have on 
EoR. 

4.2. Correlation Length Between Ionized Regions 



We measure the 3D power spectrum of the ionization 
field for each EoR model. Here the ionization field is cal- 
culated by stepping through redshift and querying each 
cell ^re(x) > z (cf. Fig. |31 bottom panels). The power 
spectrum is expected to peak on scales where the ion- 
ized regions are the most correlated. We calculate the 
wavenumber, k e , where the power spectrum peaks, and 
define the typical correlation length between ionized re- 
gions as A c ~ 2ir/k e . 

The value of A e is a proxy for the typical ionized bubble 
size. We find that the largest A c values appear around 
z. After z, A c no longer measures the typical correla- 
tion length between ionized regions, instead it measures 
the typical correlation between neutral regions (i.e. typ- 
ical size of neutral clouds). The redshift evolution of A c 
for fixed k a and a is independent of z (cf. Fig. At 
a fixed z = 10 the value for A c ranges from approxi- 
mately 3 — 90 Mpc//i between our long and short dura- 
tion reionization models (cf. Fig. [9J . Complementary to 
the results shown in Sec. 14.11 we find that the long du- 
ration reionization model has small correlations lengths 
between ionized regions, while the short duration reion- 
ization model has large correlation lengths between ion- 
ized regions. We chose to compare our models at z when 
studying the parameter space of k and a, since at z 
all these models have approximately the same ionization 
fraction (~ 50%). In FigureHUlwe show parameter space 
grid for A c a function of k and a for a fixed z = 10. 
The same trends for A z as a function of k Q and a arc 
found in A c . However, the degenerate combinations of 
k a and a that give similar A z values have different A c 
values. Thus, the model degeneracies in any observation 
that measures A z will be removed by any observation 
that measures A e such as the 21cm signal or kSZ power 
spectrum. 

5. DISCUSION 

Semi-analytic models are clearly important tools for 
understanding the EoR. They are able to quickly ex- 
plore the large parameter space of unknown and uncon- 
strained physics in an attempt to quantify how these 
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Fig. 9. — The typical correlation length between ionized re- 
gions, A e , for three fiducial parameter models with z = 8,10,12 
(green, red, and blue lines, respectively) and the long and short 
duration reionization models (orange and light green dashed line, 
respectively). The values of \ e (z) for the fiducial models is inde- 
pendent of z. Comparing the values at z = 10, where all models 
with z = 10 have similar x e values, there is a significant difference 
between X e {z) for the long and short duration reionization. 

unknown physical processes impact EoR observables. 
There is an abundan c e of semi-analyti c models (e.g. 
Furlanetto et al.1 120041 iZahn et al.1 l2007t lAlvarez et all 
2009tlChoudhurv et al1l2009h that "compute an ionization 
field from a density field (initial conditions or N-body), 
which rely on a couple of free parameters such as an ef- 
ficiency £ to model the unknown and unresolved physics 
of EoR. When comparing against simulations, these free 
parameters are tuned such that the models match the 
ionization fractions of the radiative transfers simulations. 
However, one cannot directly compare the parameter val- 
ues in the model to all values in the simulation and once 
these parameters are tuned to match simulations their 
physical interpretations diminish. Our model differs in 
that all the complex non-linear physics within the sim- 
ulations is incapsulatcd in the parametric form b mz (k). 
The direct comparison between model and simulations 
is trivial, since 6 mz (fc) is computed from the simulations 
and inserted into our model. In principle, when we cali- 
brated b mz (k) off simulations there is no loss of accuracy 
on large scales (> 1 Mpc/h) and there is only one free 
parameter, z. Additionally, our model is extremely quick 
and easily applicable to large N-body simulations, since 
it requires only two FFTs to calculate an ionization field, 
with the rate limiting step being the input and output of 
the large density and ionization fields, respectively. 

The model we proposed is based on an inside-out sce- 
nario for reionization, like several other semi-analytic 
models, which assumes that the first generation of galax- 
ies are the dominant sources for the ionizing photons. 
This model does not capture exotic reionization scenar- 
ios, like those where void regions reionize first. Fur- 
thermore, the simulations that our model is derived 
from do not include ionizing photons from population 
III s tars or first X-ray sources like high redsh ift quasars 
(e.g. IWvithe k Ce^l 120071 iTrac k CeiJl200l . The im- 
pact of these alternate sources for ion izing photons on 
the EoR is still an open question (e.g. lAhn et all 120121 : 
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Fig. 10. — Parameter grid of k and a for the typical correlation 
length between ionized regions, A c for z = 10 and at z = z. The 
previous degenerate combinations of k and a that gave the same 
values for A z give different values for A c . 

IVisbal k Loebl I2012T: IPeng et all [2011 . In future work, 
we look into implementing higher order b mz (k) and r mz 
statistics, which should increase the ability of our model 
to accurately map simulations onto larger scales. 

6. CONCLUSIONS 

We present a novel semi-analytic model for calculat- 
ing the redshift evolution of the ionization field dur- 
ing EoR that is fast and easily applicable to large vol- 
ume N-body simulations. Our model is motivated by 
and calibrated off of RadHydro simulations, which show 
there is a strong correlation between the density and the 
rcionization-rcdshift fields on scales > 1 Mpc/h. A sim- 
ple filter (cf. Eq. [7} is convolved with a non-linear den- 
sity field at z — z to obtain a zre(x), which depends 
mainly on the parameters k a , a and z. The number of 
parameters is reduced to one, z (essentially sets the mean 
rcionization-redshift), when the values of k and a are fit 
to simulation results. 

We found that this model performed well on large 
scales when we compare it directly to RadHydro simula- 
tions. Three of the comparisons we preformed were strin- 
gent tests of the simulated and modeled zre(x) fields: 
(1) We compared slices of the 2re(x) field between the 
simulation and model illustrating the minor differences 
in the evolution of the zre(x) field. (2) We compared 
the cumulative distribution of the zre(x) values, which 
again yielded minor differences between the simulation 
and model. (3) We constructed a 2D probability distri- 
bution function for a cell by cell comparison of zre(x) 
between the simulation and model which showed that 
90% of all zre(x) values in the model were within 10% 
of the simulation values. Given that all the slight differ- 
ence between the simulation and model were well below 
the uncertainty in details of the astrophysical processes 
at work during EoR, the model we proposed is a great 
tool for incorporating radiation hydrodynamic physics of 
reionization into large N-body simulations. 



10 



Battaglia et al. 



We introduce a physical understanding of the param- 
eters k a and a by comparing the ionization histories, x e , 
and typical correlation length between ionized regions, 
A e , for various combination of k a and a. Decreasing the 
parameter k Q shortens the duration of EoR, while in- 
creasing k a lengthens the duration of EoR, which is a di- 
rect result of changing the variance of S z . The parameter 
a has the opposite behavior, i.e. increasing a shortens 
the duration of EoR and vise versa. For the values of 
z we explored, these physical interpretations of k a and 
a are independent of the z value. There are degenerate 
combinations of k and a that produce nearly identical 
x c (z). These degeneracies are broken by A e and allow one 
to further differentiate between parameter combinations. 
Thus, degeneracy in observablcs, which depend on x e (z), 
can be broken when combined with obscrvables that de- 
pend on the A e . Similar to x e (z), the values for X e (z) at 
fixed k a and a arc practically independent of z. It is nec- 
essary to compare values for A (z) for varying k Q and a 
at fixed x e , so comparing at at z = z is a natural choice. 
For a fixed z = z, we demonstrate that smaller A e are 
obtained by increasing k a or decreasing a while larger 
A e are obtained by decreasing k a or increasing a, with 
values for A ranging from ~ 3 — 90 Mpc//i. In summary, 
any combination of k and a that extends the function 
bmz(k) to large values of k will increases the amount of 
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Our method is an accurate and fast tool for exploring 
galactic rcionization of large scales and going forward we 
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